Install libraries¶

In [ ]:
# from google.colab import drive
# drive.mount('/content/drive')
In [ ]:
!pip install -e .

Obtaining file:///Users/M286333/Documents/_projects/hifuku
  Preparing metadata (setup.py) ... done
Requirement already satisfied: opencv-python in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from hifuku==0.0.1) (4.6.0)
Requirement already satisfied: albumentations in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from hifuku==0.0.1) (1.3.0)
Requirement already satisfied: pytorch-lightning in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from hifuku==0.0.1) (2.0.0)
Requirement already satisfied: segmentation-models-pytorch in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from hifuku==0.0.1) (0.3.1)
Requirement already satisfied: scikit-image in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from hifuku==0.0.1) (0.19.3)
Requirement already satisfied: opencv-python-headless>=4.1.1 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from albumentations->hifuku==0.0.1) (4.7.0.72)
Requirement already satisfied: qudida>=0.0.4 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from albumentations->hifuku==0.0.1) (0.0.4)
Requirement already satisfied: scipy in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from albumentations->hifuku==0.0.1) (1.9.3)
Requirement already satisfied: numpy>=1.11.1 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from albumentations->hifuku==0.0.1) (1.23.5)
Requirement already satisfied: PyYAML in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from albumentations->hifuku==0.0.1) (6.0)
Requirement already satisfied: imageio>=2.4.1 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from scikit-image->hifuku==0.0.1) (2.23.0)
Requirement already satisfied: PyWavelets>=1.1.1 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from scikit-image->hifuku==0.0.1) (1.4.1)
Requirement already satisfied: tifffile>=2019.7.26 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from scikit-image->hifuku==0.0.1) (2022.10.10)
Requirement already satisfied: pillow!=7.1.0,!=7.1.1,!=8.3.0,>=6.1.0 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from scikit-image->hifuku==0.0.1) (9.4.0)
Requirement already satisfied: networkx>=2.2 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from scikit-image->hifuku==0.0.1) (2.8.8)
Requirement already satisfied: packaging>=20.0 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from scikit-image->hifuku==0.0.1) (22.0)
Requirement already satisfied: fsspec[http]>2021.06.0 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from pytorch-lightning->hifuku==0.0.1) (2022.11.0)
Requirement already satisfied: torchmetrics>=0.7.0 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from pytorch-lightning->hifuku==0.0.1) (0.11.0)
Requirement already satisfied: torch>=1.11.0 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from pytorch-lightning->hifuku==0.0.1) (1.13.0)
Requirement already satisfied: lightning-utilities>=0.7.0 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from pytorch-lightning->hifuku==0.0.1) (0.8.0)
Requirement already satisfied: typing-extensions>=4.0.0 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from pytorch-lightning->hifuku==0.0.1) (4.4.0)
Requirement already satisfied: tqdm>=4.57.0 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from pytorch-lightning->hifuku==0.0.1) (4.64.1)
Requirement already satisfied: timm==0.4.12 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from segmentation-models-pytorch->hifuku==0.0.1) (0.4.12)
Requirement already satisfied: torchvision>=0.5.0 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from segmentation-models-pytorch->hifuku==0.0.1) (0.14.0a0)
Requirement already satisfied: efficientnet-pytorch==0.7.1 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from segmentation-models-pytorch->hifuku==0.0.1) (0.7.1)
Requirement already satisfied: pretrainedmodels==0.7.4 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from segmentation-models-pytorch->hifuku==0.0.1) (0.7.4)
Requirement already satisfied: munch in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from pretrainedmodels==0.7.4->segmentation-models-pytorch->hifuku==0.0.1) (2.5.0)
Requirement already satisfied: requests in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from fsspec[http]>2021.06.0->pytorch-lightning->hifuku==0.0.1) (2.28.2)
Requirement already satisfied: aiohttp!=4.0.0a0,!=4.0.0a1 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from fsspec[http]>2021.06.0->pytorch-lightning->hifuku==0.0.1) (3.8.3)
Requirement already satisfied: scikit-learn>=0.19.1 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from qudida>=0.0.4->albumentations->hifuku==0.0.1) (1.2.0)
Requirement already satisfied: yarl<2.0,>=1.0 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from aiohttp!=4.0.0a0,!=4.0.0a1->fsspec[http]>2021.06.0->pytorch-lightning->hifuku==0.0.1) (1.8.2)
Requirement already satisfied: multidict<7.0,>=4.5 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from aiohttp!=4.0.0a0,!=4.0.0a1->fsspec[http]>2021.06.0->pytorch-lightning->hifuku==0.0.1) (6.0.4)
Requirement already satisfied: async-timeout<5.0,>=4.0.0a3 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from aiohttp!=4.0.0a0,!=4.0.0a1->fsspec[http]>2021.06.0->pytorch-lightning->hifuku==0.0.1) (4.0.2)
Requirement already satisfied: charset-normalizer<3.0,>=2.0 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from aiohttp!=4.0.0a0,!=4.0.0a1->fsspec[http]>2021.06.0->pytorch-lightning->hifuku==0.0.1) (2.1.1)
Requirement already satisfied: attrs>=17.3.0 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from aiohttp!=4.0.0a0,!=4.0.0a1->fsspec[http]>2021.06.0->pytorch-lightning->hifuku==0.0.1) (22.2.0)
Requirement already satisfied: frozenlist>=1.1.1 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from aiohttp!=4.0.0a0,!=4.0.0a1->fsspec[http]>2021.06.0->pytorch-lightning->hifuku==0.0.1) (1.3.3)
Requirement already satisfied: aiosignal>=1.1.2 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from aiohttp!=4.0.0a0,!=4.0.0a1->fsspec[http]>2021.06.0->pytorch-lightning->hifuku==0.0.1) (1.3.1)
Requirement already satisfied: threadpoolctl>=2.0.0 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from scikit-learn>=0.19.1->qudida>=0.0.4->albumentations->hifuku==0.0.1) (3.1.0)
Requirement already satisfied: joblib>=1.1.1 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from scikit-learn>=0.19.1->qudida>=0.0.4->albumentations->hifuku==0.0.1) (1.2.0)
Requirement already satisfied: six in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from munch->pretrainedmodels==0.7.4->segmentation-models-pytorch->hifuku==0.0.1) (1.16.0)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from requests->fsspec[http]>2021.06.0->pytorch-lightning->hifuku==0.0.1) (1.26.13)
Requirement already satisfied: idna<4,>=2.5 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from requests->fsspec[http]>2021.06.0->pytorch-lightning->hifuku==0.0.1) (3.4)
Requirement already satisfied: certifi>=2017.4.17 in /Users/M286333/miniforge3/envs/py39/lib/python3.9/site-packages (from requests->fsspec[http]>2021.06.0->pytorch-lightning->hifuku==0.0.1) (2023.5.7)
Installing collected packages: hifuku
  Attempting uninstall: hifuku
    Found existing installation: hifuku 0.0.1
    Uninstalling hifuku-0.0.1:
      Successfully uninstalled hifuku-0.0.1
  Running setup.py develop for hifuku
Successfully installed hifuku-0.0.1
In [ ]:
# %cd /content/drive/MyDrive/hifuku
# !pip install -r requirements.txt
In [ ]:
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

import os
from glob import glob

import csv
import os.path
from pathlib import Path
import random
import math
import warnings
warnings.simplefilter('ignore')

import cv2
import PIL
PIL.Image.MAX_IMAGE_PIXELS = 933120000
from PIL import Image

import joblib
from sklearn.manifold import TSNE
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler
from sklearn.mixture import GaussianMixture
from sklearn.cluster import KMeans
from scipy import stats
import ast
from scipy.spatial import distance
# from scipy.spatial import Voronoi, voronoi_plot_2d
from scipy.stats import entropy
from sklearn.neighbors import KernelDensity

from shapely.geometry import Polygon
from tqdm import tqdm

import hifuku
In [ ]:
palette = sns.color_palette('viridis', n_colors=5)
sns.set(context='talk', style='ticks', palette=palette, font_scale=.7,rc={
    'axes.spines.right': False,
    'axes.spines.top': False,
    'figure.figsize': [4.0, 4.0],
    'legend.fancybox': False,
    'legend.edgecolor': 'black',
    'legend.fancybox': False,
    'legend.frameon': True,
    'legend.handletextpad': 0.5,
    'legend.loc': 'upper right',
    'legend.loc': 'best',
    'legend.title_fontsize': None,
})

sns.palplot(palette)
colors = [(int(color[0]*255), int(color[1]*255), int(color[2]*255)) for color in palette]

Set parameters¶

In [ ]:
# root = '/content/drive/MyDrive/hifuku'
root = '/Users/M286333/Documents/_projects/hifuku'

# a sample WSI image
is_wsi = True
path = f'{root}/data/sample_wsi_scale_0273.jpg'
scale = 0.2730 # µm / pixel


# if you want to analyze a non-wsi sample image, activate the following 3 lines 
# is_wsi = False
# path = f'{root}/data/sample_scale_02325.jpg'
# scale = 0.2325 # µm / pixel


image_id = os.path.splitext(os.path.splitext(os.path.basename(path))[0])[0]
save_dir= f'{root}/results/{image_id}'
image_id
Out[ ]:
'sample_wsi_scale_0273'

Run Hifuku¶

In [ ]:
hifuku.main(root, path, is_wsi=is_wsi, scale=scale)

Global seed set to 0

<Figure size 400x400 with 0 Axes>
<Figure size 2000x2000 with 0 Axes>

See results¶

In [ ]:
df_nerve = pd.read_csv(f'{root}/results/{image_id}/data_nerve.csv', index_col=0)
df_nerve
Out[ ]:
id is_wsi x_px y_px scale total_fas n_fibers total_area total_density std_density
0 sample_wsi_scale_0273 True 7813 7863 0.273 10 4155 0.583923 7115.665869 1005.693849
In [ ]:
df_fas = pd.read_csv(f'{root}/results/{image_id}/data_fas.csv', index_col=0)
df_fas
Out[ ]:
n_fas n_fib area density
0 0 1119 0.182891 6118.394700
1 1 383 0.050129 7640.306931
2 2 946 0.140433 6736.305686
3 3 534 0.065160 8195.160650
4 4 366 0.043105 8490.945248
5 5 125 0.015002 8332.047390
6 6 164 0.023515 6974.354989
7 7 337 0.036705 9181.195894
8 8 106 0.016642 6369.273575
9 9 75 0.010340 7253.538672
In [ ]:
df_fib = pd.read_csv(f'{root}/results/{image_id}/data_fib.csv', index_col=0)
df_fib = df_fib.rename(columns={'diameter_out': 'diameter'})
df_fib
Out[ ]:
area_out area_in perimeter circularity convexity solidity eccentricity major_axis_length minor_axis_length angle diameter diameter_in thickness g_ratio n_fas xmin ymin xmax ymax
0 7.677978 3.292691 231.036577 0.606331 0.933163 0.906069 2.445910 5.186807 2.120604 34.596996 3.126644 2.047531 0.539556 0.654865 0.0 880.977778 205.297297 901.511111 227.675676
1 118.794754 39.889411 830.482315 0.726041 0.897158 0.973624 1.431737 14.793043 10.332236 97.859131 12.298544 7.126624 2.585960 0.579469 0.0 892.711111 172.216216 953.333333 218.918919
2 18.742553 9.934716 370.676186 0.574995 0.876249 0.875871 1.958578 7.144187 3.647640 174.681580 4.885055 3.556582 0.664237 0.728054 0.0 983.644444 191.675676 1006.133333 223.783784
3 105.212589 33.983733 748.825463 0.790917 0.916607 0.979055 1.222804 12.935622 10.578653 123.923126 11.574145 6.577951 2.498097 0.568331 0.0 929.866667 141.081081 980.711111 188.756757
4 27.972224 15.121934 386.475176 0.789420 0.921415 0.967320 1.252309 6.724402 5.369602 130.612839 5.967859 4.387920 0.789969 0.735259 0.0 1129.333333 154.702703 1159.644444 184.864865
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
4150 32.011696 8.365135 413.161468 0.790485 0.932958 0.977115 1.480821 7.836022 5.291674 131.498795 6.384243 3.263560 1.560342 0.511190 9.0 217.846154 320.307692 250.153846 351.692308
4151 46.576153 14.500362 502.173661 0.778539 0.932473 0.970675 1.495623 9.482864 6.340413 105.341042 7.700818 4.296794 1.702012 0.557966 9.0 422.769231 267.692308 463.384615 298.153846
4152 47.135121 17.149123 508.960457 0.767010 0.927581 0.956706 1.400718 9.372693 6.691348 97.500267 7.746890 4.672787 1.537051 0.603182 9.0 424.615385 237.230769 463.384615 268.615385
4153 4.652100 1.997377 169.095453 0.685819 0.920109 0.919022 1.585113 3.185095 2.009380 93.401825 2.433770 1.594722 0.419524 0.655248 9.0 430.153846 230.769231 448.615385 244.615385
4154 9.256502 4.119963 262.509666 0.566214 0.899516 0.936369 2.530454 5.572750 2.202273 88.055771 3.433037 2.290349 0.571344 0.667150 9.0 257.538462 390.461538 283.384615 408.000000

4155 rows × 19 columns

In [ ]:
img = Image.open(f'{save_dir}/masked_fascicles.jpg')
plt.imshow(img, aspect='auto')
plt.axis('off')
Out[ ]:
(-0.5, 7862.5, 7812.5, -0.5)
In [ ]:
col = 3
row = len(df_fas)//col + 1
plt.figure(figsize=(col*3, row*3))
plt.subplots_adjust(wspace=0.1, hspace=0.1)

for num in range(len(df_fas)):
    plt.subplot(row, col, num+1)
    img = Image.open(f'{save_dir}/fascicles/fascicle_{num:03}_bbox.jpg')
    plt.imshow(img)
    # plt.title(f'fascicle {num}', fontsize=10)
    plt.axis('off')

Exclude improperly detected fascicle(s), if present.¶

In [ ]:
fas_to_remove = [] # fill number(s) of improperly detected fas in the list
df_fib_rm = df_fib[~df_fib['n_fas'].isin(fas_to_remove)].dropna(axis=0)
df_fib_rm
Out[ ]:
area_out area_in perimeter circularity convexity solidity eccentricity major_axis_length minor_axis_length angle diameter diameter_in thickness g_ratio n_fas xmin ymin xmax ymax
0 7.677978 3.292691 231.036577 0.606331 0.933163 0.906069 2.445910 5.186807 2.120604 34.596996 3.126644 2.047531 0.539556 0.654865 0.0 880.977778 205.297297 901.511111 227.675676
1 118.794754 39.889411 830.482315 0.726041 0.897158 0.973624 1.431737 14.793043 10.332236 97.859131 12.298544 7.126624 2.585960 0.579469 0.0 892.711111 172.216216 953.333333 218.918919
2 18.742553 9.934716 370.676186 0.574995 0.876249 0.875871 1.958578 7.144187 3.647640 174.681580 4.885055 3.556582 0.664237 0.728054 0.0 983.644444 191.675676 1006.133333 223.783784
3 105.212589 33.983733 748.825463 0.790917 0.916607 0.979055 1.222804 12.935622 10.578653 123.923126 11.574145 6.577951 2.498097 0.568331 0.0 929.866667 141.081081 980.711111 188.756757
4 27.972224 15.121934 386.475176 0.789420 0.921415 0.967320 1.252309 6.724402 5.369602 130.612839 5.967859 4.387920 0.789969 0.735259 0.0 1129.333333 154.702703 1159.644444 184.864865
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
4150 32.011696 8.365135 413.161468 0.790485 0.932958 0.977115 1.480821 7.836022 5.291674 131.498795 6.384243 3.263560 1.560342 0.511190 9.0 217.846154 320.307692 250.153846 351.692308
4151 46.576153 14.500362 502.173661 0.778539 0.932473 0.970675 1.495623 9.482864 6.340413 105.341042 7.700818 4.296794 1.702012 0.557966 9.0 422.769231 267.692308 463.384615 298.153846
4152 47.135121 17.149123 508.960457 0.767010 0.927581 0.956706 1.400718 9.372693 6.691348 97.500267 7.746890 4.672787 1.537051 0.603182 9.0 424.615385 237.230769 463.384615 268.615385
4153 4.652100 1.997377 169.095453 0.685819 0.920109 0.919022 1.585113 3.185095 2.009380 93.401825 2.433770 1.594722 0.419524 0.655248 9.0 430.153846 230.769231 448.615385 244.615385
4154 9.256502 4.119963 262.509666 0.566214 0.899516 0.936369 2.530454 5.572750 2.202273 88.055771 3.433037 2.290349 0.571344 0.667150 9.0 257.538462 390.461538 283.384615 408.000000

3893 rows × 19 columns

In [ ]:
len(df_fib_rm[df_fib_rm['diameter'] >= 25]), df_fib_rm.isna().sum().sum()
Out[ ]:
(0, 0)
In [ ]:
df_fas_rm = df_fas[~df_fas['n_fas'].isin(fas_to_remove)]
df_fas_rm
Out[ ]:
n_fas n_fib area density
0 0 1119 0.182891 6118.394700
1 1 383 0.050129 7640.306931
2 2 946 0.140433 6736.305686
3 3 534 0.065160 8195.160650
4 4 366 0.043105 8490.945248
5 5 125 0.015002 8332.047390
6 6 164 0.023515 6974.354989
7 7 337 0.036705 9181.195894
8 8 106 0.016642 6369.273575
9 9 75 0.010340 7253.538672
In [ ]:
df_nerve_rm = df_nerve.copy()
df_nerve_rm['total_fas'] = len(df_fas_rm)
df_nerve_rm['n_fibers'] = df_fas_rm['n_fib'].sum()
df_nerve_rm['total_area'] = df_fas_rm['area'].sum()
df_nerve_rm['total_density'] = df_nerve_rm['n_fibers'] / df_nerve_rm['total_area']
df_nerve_rm['std_density'] = df_fas_rm['density'].std()
df_nerve_rm
Out[ ]:
id is_wsi x_px y_px scale total_fas n_fibers total_area total_density std_density
0 sample_wsi_scale_0273 True 7813 7863 0.273 10 4155 0.583923 7115.665869 1005.693849

Visualize results¶

In [ ]:
def get_density(df, area, name, xlim=20, dev=1):
    idx = []
    density = []
    bins = int(xlim // dev) + 1
    for i in range(bins):
        if i == bins - 1:
            dens = len(df[df[name] >= (i) * dev]) / area
        else:
            dens = len(df[ (df[name] >= (i) * dev) & (df[name] < (i + 1) * dev) ]) / area
        idx.append(i * dev)
        density.append(dens)
    density_df = pd.DataFrame({name: density}, index = idx)
    return density_df

def draw_histogram(df, total_area, name='diameter', dev=1.0, width=1, xlim=20, tick=1):
    density = get_density(df, total_area, name=name, xlim=xlim, dev=dev)
    fig, ax = plt.subplots(figsize=(6, 3))

    # ax.set_title(name)
    ax.set_xlim(0, xlim)
    ax.set_xlabel(f'{name} [um]', fontsize=10)
    ax.set_ylabel('Density [/mm2]', fontsize=10)
    ax.set_xticks(np.arange(0, xlim, tick))
    ax.set_xticks(np.arange(0, xlim, dev), minor=True)


    ax.bar(density.index, density[name], width=width, color='teal')
    # ax.bar(density_ref.index+width/2, density_ref[name], width=width, alpha=0.5, color='powderblue')
    # ax.legend(['sample', 'normal 57 y.o. male'])
    plt.savefig(f'{save_dir}/histogram_{name}.jpg')
In [ ]:
area = df_nerve_rm['total_area'][0]
draw_histogram(df_fib_rm, area, name='diameter', dev=1, width=1, xlim=20, tick=1)
draw_histogram(df_fib_rm, area, name='thickness', dev=0.20, width=0.20, xlim=5, tick=1)
In [ ]:
palette = sns.color_palette('gray', n_colors=2)
sns.set(context='talk', style='ticks', palette=palette, font_scale=.7,rc={
    'axes.spines.right': False,
    'axes.spines.top': False,
    'figure.figsize': [4.0, 4.0],
    'legend.fancybox': False,
    'legend.edgecolor': 'black',
    'legend.fancybox': False,
    'legend.frameon': True,
    'legend.handletextpad': 0.5,
    'legend.loc': 'upper right',
    'legend.loc': 'best',
    'legend.title_fontsize': None,
})
sns.palplot(palette)
In [ ]:
sns.histplot(df_fib_rm['diameter'], bins=40)
plt.xlim(0, 20)
plt.xlabel('Diameter [um]', fontsize=10)
plt.ylabel('Count / slide', fontsize=10)
plt.savefig(f'{save_dir}/histogram_diameter.jpg')
In [ ]:
sns.histplot(df_fib_rm['thickness'], bins=50)
plt.xlim(0, 5)
plt.xlabel('Thickness [um]', fontsize=10)
plt.ylabel('Count / slide', fontsize=10)
plt.savefig(f'{save_dir}/histogram_thickness.jpg')
In [ ]:
sns.scatterplot(data=df_fib_rm, x='diameter', y='thickness', s=2)
plt.xlabel('Diameter [um]', fontsize=10)
plt.ylabel('Thickness [um]', fontsize=10)
plt.xlim(0, 20)
plt.ylim(0, 5)
plt.savefig(f'{save_dir}/scatter_diameter_thickness.jpg')
In [ ]:
sns.scatterplot(data=df_fib_rm, x='diameter', y='g_ratio', s=2)
plt.xlim(0, 20)
plt.ylim(0, 1)
plt.savefig(f'{save_dir}/scatter_diameter_g_ratio.jpg')

Classify large and small fibers with Gaussian Mixture model¶

In [ ]:
gmm = GaussianMixture(n_components=2, covariance_type='tied')
loaded_gmm = joblib.load(f'{root}/weights/gmm_classifier.pkl')
df_fib_rm['gmm'] = loaded_gmm.predict(df_fib_rm[['diameter', 'thickness']])
In [ ]:
palette = sns.color_palette('viridis', n_colors=3)
sns.set(context='talk', style='ticks', palette=palette, font_scale=.7,rc={
    'axes.spines.right': False,
    'axes.spines.top': False,
    'figure.figsize': [4.0, 4.0],
    'legend.fancybox': False,
    'legend.edgecolor': 'black',
    'legend.fancybox': False,
    'legend.frameon': True,
    'legend.handletextpad': 0.5,
    'legend.loc': 'upper right',
    'legend.loc': 'best',
    'legend.title_fontsize': None,
})

sns.palplot(palette)
In [ ]:
sns.scatterplot(data=df_fib_rm, x='diameter', y='g_ratio', hue='gmm', s=2)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.xlabel('Diameter [um]', fontsize=10)
plt.ylabel('g-ratio', fontsize=10)
plt.xlim(0, 20)
plt.ylim(0, 1)
plt.savefig(f'{save_dir}/scatter_diameter_g_ratio_gmm.jpg')
In [ ]:
sns.scatterplot(data=df_fib_rm, x='diameter', y='thickness', hue='gmm', s=2)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.xlabel('Diameter [um]', fontsize=10)
plt.ylabel('Thickness [um]', fontsize=10)
plt.xlim(0, 20)
plt.ylim(0, 5)
plt.savefig(f'{save_dir}/scatter_diameter_thickness_gmm.jpg')
In [ ]:
sns.histplot(data=df_fib_rm, x='diameter', bins=40, hue='gmm')
labels = ['small fiber', 'large fiber']
plt.legend(labels=labels, bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.xlim(0, 20)
plt.xlabel('Diameter [um]', fontsize=10)
plt.ylabel('Count / slide', fontsize=10)
plt.savefig(f'{save_dir}/histogram_diameter_gmm.jpg')
In [ ]:
sns.histplot(data=df_fib_rm, x='thickness', bins=50, hue='gmm')
labels = ['small fiber', 'large fiber']
plt.legend(labels=labels, bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.xlim(0, 5)
plt.xlabel('Thickness [um]', fontsize=10)
plt.ylabel('Count / slide', fontsize=10)
plt.savefig(f'{save_dir}/histogram_thickness_gmm.jpg')

Calculate linear regression¶

In [ ]:
def get_linear_regression(df, xlabel='diameter', ylabel='g_ratio'):
    df = df.copy().dropna(subset=(xlabel, ylabel))
    x = df[xlabel].values.reshape(-1, 1)
    y = df[ylabel]
    lr = LinearRegression()
    lr.fit(x, y)
    # r2 = lr.score(x, y)
    coef = lr.coef_
    intercept = lr.intercept_


    sns.jointplot(data=df_fib_rm, x='diameter', xlim=(0, 15), ylim=(0, 1), y='g_ratio', hue='gmm', s=2)
    plt.legend([],[], frameon=False)

    x_ = np.arange(min(df[xlabel]), max(df[xlabel]), 0.1)
    plt.plot(x_, coef*x_+intercept, color='gray', linewidth=2)



    plt.xlim(0, 20)
    plt.ylim(0, 1)
    plt.xlabel('Diameter [um]', fontsize=10)
    plt.ylabel('g-ratio', fontsize=10)

    plt.savefig(f'{save_dir}/{xlabel}_{ylabel}_lr.jpg')
    print(
        f'coef:{coef[0]:.4f}, intercept:{intercept:.3f}'
    )
    return coef[0], intercept
In [ ]:
df_nerve['coef'], df_nerve['intercept'] = get_linear_regression(df_fib_rm)

coef:-0.0156, intercept:0.723
In [ ]:
df_fib
Out[ ]:
area_out area_in perimeter circularity convexity solidity eccentricity major_axis_length minor_axis_length angle diameter diameter_in thickness g_ratio n_fas xmin ymin xmax ymax
0 7.677978 3.292691 231.036577 0.606331 0.933163 0.906069 2.445910 5.186807 2.120604 34.596996 3.126644 2.047531 0.539556 0.654865 0.0 880.977778 205.297297 901.511111 227.675676
1 118.794754 39.889411 830.482315 0.726041 0.897158 0.973624 1.431737 14.793043 10.332236 97.859131 12.298544 7.126624 2.585960 0.579469 0.0 892.711111 172.216216 953.333333 218.918919
2 18.742553 9.934716 370.676186 0.574995 0.876249 0.875871 1.958578 7.144187 3.647640 174.681580 4.885055 3.556582 0.664237 0.728054 0.0 983.644444 191.675676 1006.133333 223.783784
3 105.212589 33.983733 748.825463 0.790917 0.916607 0.979055 1.222804 12.935622 10.578653 123.923126 11.574145 6.577951 2.498097 0.568331 0.0 929.866667 141.081081 980.711111 188.756757
4 27.972224 15.121934 386.475176 0.789420 0.921415 0.967320 1.252309 6.724402 5.369602 130.612839 5.967859 4.387920 0.789969 0.735259 0.0 1129.333333 154.702703 1159.644444 184.864865
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
4150 32.011696 8.365135 413.161468 0.790485 0.932958 0.977115 1.480821 7.836022 5.291674 131.498795 6.384243 3.263560 1.560342 0.511190 9.0 217.846154 320.307692 250.153846 351.692308
4151 46.576153 14.500362 502.173661 0.778539 0.932473 0.970675 1.495623 9.482864 6.340413 105.341042 7.700818 4.296794 1.702012 0.557966 9.0 422.769231 267.692308 463.384615 298.153846
4152 47.135121 17.149123 508.960457 0.767010 0.927581 0.956706 1.400718 9.372693 6.691348 97.500267 7.746890 4.672787 1.537051 0.603182 9.0 424.615385 237.230769 463.384615 268.615385
4153 4.652100 1.997377 169.095453 0.685819 0.920109 0.919022 1.585113 3.185095 2.009380 93.401825 2.433770 1.594722 0.419524 0.655248 9.0 430.153846 230.769231 448.615385 244.615385
4154 9.256502 4.119963 262.509666 0.566214 0.899516 0.936369 2.530454 5.572750 2.202273 88.055771 3.433037 2.290349 0.571344 0.667150 9.0 257.538462 390.461538 283.384615 408.000000

4155 rows × 19 columns

Spatial analysis¶

In [ ]:
# down scale to 1/10 to compute fast
df_fib['x'] = (df_fib['xmin'] + df_fib['xmax']) / 20
df_fib['y'] = (df_fib['ymin'] + df_fib['ymax']) / 20

df_fas_rm['NND'] = np.nan
df_fas_rm['KDE_entropy'] = np.nan
for idx in range(len(df_fas_rm)):
    data = df_fib[(df_fib['n_fas']==idx)].dropna(subset=['x', 'y'])
    data = data[['x', 'y']]
    x = data['x']
    y = data['y']
    if len(data) < 10:
        continue

    # Calculate nearest neighbor distances
    distances = distance.cdist(data, data)
    nnd = np.min(np.ma.masked_array(distances, mask=np.eye(len(data))), axis=1)
    # Find the minimum distance for each point
    masked_distances = np.ma.masked_array(distances, mask=np.eye(len(data)))
    nearest_distances = masked_distances.min(axis=1)
    # Calculate the mean nearest neighbor distance
    mean_nnd = nearest_distances.mean() * 10 * scale # um
    df_fas_rm['NND'][idx] =  mean_nnd
    print(f"Mean Nearest Neighbor Distance: {mean_nnd:.2f} um")


    # Perform kernel density estimation
    bandwidth = 5  # note that 1/10 scale to original pixel
    kde = KernelDensity(bandwidth=bandwidth, kernel='gaussian')
    kde.fit(data)

    # Evaluate the KDE on the grid
    grid_size = 1
    grid_x, grid_y = np.mgrid[0:x.max():grid_size, 0:y.max():grid_size]
    grid_points = np.column_stack((grid_x.ravel(), grid_y.ravel()))
    log_density = kde.score_samples(grid_points)
    density = np.exp(log_density)

    # Calculate the entropy of the KDE distribution
    entropy_value = entropy(density)
    df_fas_rm['KDE_entropy'][idx] =  entropy_value
    print(f"Entropy: {entropy_value:.2f}")

    # plot data
    fig, ax = plt.subplots(1, 4, figsize=(12, 3))
    img = Image.open(f'{save_dir}/fascicles/fascicle_{idx:03}_bbox.jpg')
    ax[0].imshow(img)
    ax[1].set_aspect('equal')

    ax[1].scatter(x, y, s=5)
    ax[1].set_xlim(0, x.max()+10)
    ax[1].set_ylim(0, y.max()+10)
    ax[1].invert_yaxis()
    ax[1].axis('off')
    ax[1].set_aspect('equal')

    im = ax[2].imshow(density.reshape(grid_x.shape).T, origin='lower', cmap='viridis')
    ax[2].invert_yaxis()
    ax[2].axis('off')
    ax[2].set_aspect('equal')

    # colorbar
    cbar = fig.colorbar(im, ax=ax[3], pad=0.05, location='left')
    cbar.set_ticks([])
    ax[3].axis('off')
    plt.savefig(f'{save_dir}/fascicles/fascicle_{idx:03}_kde.jpg')
    # plt.close()

df_fas_rm['expected_nnd'] = 1 / (2 * np.sqrt(df_fas_rm['n_fib'] / (df_fas_rm['area'] * 1000000))) # um
# Calculate the ratio of observed mean nearest neighbor distance to expected value
df_fas_rm['NNI'] = df_fas_rm['NND'] / df_fas_rm['expected_nnd']

Mean Nearest Neighbor Distance: 8.13 um
Entropy: 10.10
Mean Nearest Neighbor Distance: 7.50 um
Entropy: 8.89
Mean Nearest Neighbor Distance: 8.12 um
Entropy: 9.87
Mean Nearest Neighbor Distance: 7.04 um
Entropy: 9.13
Mean Nearest Neighbor Distance: 7.16 um
Entropy: 8.75
Mean Nearest Neighbor Distance: 6.81 um
Entropy: 7.64
Mean Nearest Neighbor Distance: 7.88 um
Entropy: 8.12
Mean Nearest Neighbor Distance: 6.77 um
Entropy: 8.59
Mean Nearest Neighbor Distance: 8.01 um
Entropy: 7.86
Mean Nearest Neighbor Distance: 7.07 um
Entropy: 7.38
In [ ]:
idx = 0
# for idx in range(len(df_fas_rm)):
data = df_fib[(df_fib['n_fas']==idx)].dropna(subset=['x', 'y'])
data = data[['x', 'y']]
x = data['x']
y = data['y']
# if len(data) < 10:
#     continue

# Calculate nearest neighbor distances
distances = distance.cdist(data, data)
nnd = np.min(np.ma.masked_array(distances, mask=np.eye(len(data))), axis=1)
# Find the minimum distance for each point
masked_distances = np.ma.masked_array(distances, mask=np.eye(len(data)))
nearest_distances = masked_distances.min(axis=1)
# Calculate the mean nearest neighbor distance
mean_nnd = nearest_distances.mean() * 10 * scale # um
df_fas_rm['NND'][idx] =  mean_nnd
print(f"Mean Nearest Neighbor Distance: {mean_nnd:.2f} um")


# Perform kernel density estimation
bandwidth = 5  # note that 1/10 scale to original pixel
kde = KernelDensity(bandwidth=bandwidth, kernel='gaussian')
kde.fit(data)

# Evaluate the KDE on the grid
grid_size = 1
grid_x, grid_y = np.mgrid[0:x.max():grid_size, 0:y.max():grid_size]
grid_points = np.column_stack((grid_x.ravel(), grid_y.ravel()))
log_density = kde.score_samples(grid_points)
density = np.exp(log_density)

# Calculate the entropy of the KDE distribution
entropy_value = entropy(density)
df_fas_rm['KDE_entropy'][idx] =  entropy_value
print(f"Entropy: {entropy_value:.2f}")

# plot data
fig, ax = plt.subplots(1, 4, figsize=(12, 3))
img = Image.open(f'{save_dir}/fascicles/fascicle_{idx:03}_bbox.jpg')
ax[0].imshow(img)
ax[1].set_aspect('equal')

ax[1].scatter(x, y, s=5)
ax[1].set_xlim(0, x.max()+10)
ax[1].set_ylim(0, y.max()+10)
ax[1].invert_yaxis()
ax[1].axis('off')
ax[1].set_aspect('equal')

im = ax[2].imshow(density.reshape(grid_x.shape).T, origin='lower', cmap='viridis')
ax[2].invert_yaxis()
ax[2].axis('off')
ax[2].set_aspect('equal')

# colorbar
cbar = fig.colorbar(im, ax=ax[3], pad=0.05, location='left')
cbar.set_ticks([])
ax[3].axis('off')
plt.savefig(f'{save_dir}/fascicles/fascicle_{idx:03}_kde.jpg')
# plt.close()

Mean Nearest Neighbor Distance: 8.13 um
Entropy: 10.10
In [ ]:
density.min(), density.max(), density.mean()
Out[ ]:
(4.901244531291384e-31, 7.50410256471614e-05, 2.9131704750970153e-05)
In [ ]:
density.sum()
Out[ ]:
0.987790333117873
In [ ]:
entropy(density[:4])
Out[ ]:
0.5612423377111373
In [ ]:
density.shape
Out[ ]:
(1800,)
In [ ]:
df_fas_rm
Out[ ]:
n_fas n_fib area density NND KDE_entropy expected_nnd NNI
0 0 1119 0.182891 6118.394700 8.127743 10.104460 6.392213 1.271507
1 1 383 0.050129 7640.306931 7.495048 8.894850 5.720245 1.310267
2 2 946 0.140433 6736.305686 8.118978 9.869836 6.091989 1.332730
3 3 534 0.065160 8195.160650 7.043000 9.134708 5.523206 1.275165
4 4 366 0.043105 8490.945248 7.158880 8.746445 5.426152 1.319329
5 5 125 0.015002 8332.047390 6.814293 7.636691 5.477648 1.244018
6 6 164 0.023515 6974.354989 7.878395 8.122971 5.987120 1.315891
7 7 337 0.036705 9181.195894 6.773026 8.586368 5.218196 1.297963
8 8 106 0.016642 6369.273575 8.008600 7.862696 6.265057 1.278296
9 9 75 0.010340 7253.538672 7.070057 7.382267 5.870770 1.204281
In [ ]:
df_nerve_rm
Out[ ]:
id is_wsi x_px y_px scale total_fas n_fibers total_area total_density std_density
0 sample_wsi_scale_0273 True 7813 7863 0.273 10 4155 0.583923 7115.665869 1005.693849
In [ ]:
df_nerve_rm.to_csv(f'{save_dir}/df_nerve_rm.csv')
df_fas_rm.to_csv(f'{save_dir}/df_fas_rm.csv')
df_fib_rm.to_csv(f'{save_dir}/df_fib_rm.csv')
In [ ]:
# from google.colab import runtime
# runtime.unassign()